Last updated on 2025-09-05.
Show code
library (dplyr)library (dbplyr)library (ggplot2)for (i in list.files (here:: here ("R" ), full.names = TRUE )) {source (i)# SAMPLE NAME ## specify sample name <- c (# dmel # "dmel-head", # "dmel-body", # "dmel-heart" # # mmus "mmus-brain_stem" , "mmus-cerebellum" , "mmus-hypothalamus" , "mmus-heart" # # panu # "panu-brain_stem", # "panu-cerebellum", # "panu-hypothalamus", # "panu-heart", # "panu-scn" # sample.cycles <- c("LD", "DD") ## SPECIFY THE DATASET TO BUILD GCN WITH <- sample.names[1 ]writeLines (:: glue ("Reference tissue: {which.sample}" )
Reference tissue: mmus-brain_stem
Structure of reference tissue data:
ONE-SHOT (make_modules)
Make and plot modules
Show code
# mods <- purrr::map( # sample.names, # ~ timecourseRnaseq::make_modules( # data = tidydata.db[[.x]], # log2 = TRUE, # id_column = "gene_name", # min_expression = NULL, # automatically estimated # min_timepoints = NULL, # automatically estimated # method = "wgcna", # qc = TRUE, # sim_method = "kendall", # soft_power = 15, # automatically estimated # min_module_size = 50, # max_modules = 16, # merge_cutoff_similarity = 0.9, # plot_network_min_edge = 0.5, # plot_network = FALSE, # tidy_modules = TRUE # ) # ) |> # setNames( # sample.names # ) <- timecourseRnaseq:: make_modules (data = tidydata.db[[which.sample]],log2 = TRUE ,id_column = "gene_name" ,min_expression = NULL , # automatically estimated min_timepoints = NULL , # automatically estimated method = "wgcna" ,qc = TRUE ,sim_method = "kendall" ,soft_power = NULL , # automatically estimated min_module_size = 50 ,max_modules = 16 ,merge_cutoff_similarity = 0.9 ,plot_network = FALSE ,# plot_network_min_edge = 0.5, # only used when `plot_network == TRUE` tidy_modules = TRUE
---------------------------------------------------
1. Log2-transform and subset
---------------------------------------------------
Applying log2-transformation...Done.
Estimated min_expression = 8
Estimated min_timepoints = 6
Subsetting data...Done.
[ NOTE ]: After subsetting, 8893 of 17406 rows remain.
Flagging genes and samples with too many missing values...
..step 1
All okay!Visualizing the log-transformed data
---------------------------------------------------
2. Calculate similarity of expression
---------------------------------------------------
Running gene-gene similarity...Done.
---------------------------------------------------
3. Create adjacency matrix
---------------------------------------------------
Performing network topology analysis to pick
soft-thresholding power...
Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
1 1 0.7240 1.9500 0.654 4510 4900.0 6020
2 2 0.5090 0.5730 0.389 2970 3230.0 4790
3 3 0.0298 0.0544 0.363 2170 2280.0 4040
4 4 0.2950 -0.1720 0.710 1690 1680.0 3510
5 5 0.6060 -0.3070 0.785 1360 1280.0 3110
6 6 0.7420 -0.3970 0.817 1130 995.0 2800
7 7 0.8110 -0.4740 0.818 951 790.0 2540
8 8 0.8540 -0.5280 0.841 815 634.0 2330
9 9 0.8860 -0.5790 0.866 708 516.0 2140
10 10 0.8960 -0.6210 0.868 621 424.0 1980
11 12 0.9040 -0.6890 0.877 489 294.0 1730
12 15 0.9040 -0.7770 0.886 359 182.0 1440
13 18 0.8990 -0.8420 0.895 275 118.0 1220
14 21 0.8930 -0.8910 0.899 216 80.2 1050
Plotting resutls from the network topology analysis...
Done.
[ NOTE, FIGURE ]: Red horizontal line indicates a signed R^2 of 0.9
Setting soft-thresholding power to: 8.
Power-transforming the gene-gene similarity matrix...Done.
---------------------------------------------------
4. Convert into topological overlap matrix (dissTOM)
---------------------------------------------------
Creating dissTOM...Done.
Performing hierarchical clustering on dissTOM...Done.
---------------------------------------------------
5. Identify modules (clusters)
---------------------------------------------------
Merging modules that have a correlation ≥ 0.9 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.
Module (cluster) size:
mergedColors
antiquewhite4 bisque4 blue coral1 coral2
215 669 1938 60 287
cyan darkmagenta darkolivegreen darkorange2 darkred
337 117 435 96 157
darkseagreen4 darkslateblue floralwhite grey60 honeydew1
61 87 97 172 388
ivory lavenderblush3 lightcyan lightcyan1 lightgreen
98 155 174 98 165
magenta mediumorchid mediumpurple3 navajowhite2 orange
202 55 105 76 154
orangered4 paleturquoise palevioletred3 pink plum1
108 123 301 623 110
purple red salmon skyblue skyblue2
199 234 186 127 54
steelblue white yellow4 yellowgreen
123 146 50 111
Merging modules that have a correlation ≥ 0.85 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.
Module (cluster) size:
mergedColors
bisque4 coral1 darkgrey darkmagenta darkorange2
1397 440 3195 117 194
darkred darkseagreen4 grey60 ivory magenta
157 61 172 672 325
mediumorchid mediumpurple3 navajowhite2 orange palevioletred3
163 105 76 154 736
plum1 saddlebrown skyblue skyblue2 white
110 287 282 54 146
yellow4
50
Merging modules that have a correlation ≥ 0.8 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.
Module (cluster) size:
mergedColors
bisque4 blue coral1 darkmagenta darkorange2
2462 3195 440 796 194
darkseagreen4 ivory magenta mediumorchid mediumpurple3
61 748 325 163 105
orange skyblue2 white yellow4
154 54 146 50
Cutoff used: 0.8
Number of modules identified: 14
Calculating module-module similarity based
on module-eigengene-expression...Done.
Tidying module names...Done.
Plotting adjacency matrix for module-module similarity...
---------------------------------------------------
6. Tidy modules (clusters)
---------------------------------------------------
Modify network plot
Internal function; use ::: to call
Show code
# trash <- purrr::map( # which.sample, # ~ timecourseRnaseq:::plot_adj_as_network( # matrix = mods[[.x]][["adj_matrix_ME"]][["ME"]], # # layout = igraph::layout.sugiyama, # layout = igraph::layout_in_circle, # changed # min_edge = 0.6, # node_label_size = 1.2, # node_size = 20, # edge_size = 3, # node_frame_col = "grey20", # node_fill_col = "grey80", # vertex.frame.width = 3 # ) # ) ::: plot_adj_as_network (matrix = mods[["adj_matrix_ME" ]][["ME" ]],# layout = igraph::layout.sugiyama, layout = igraph:: layout_in_circle, # changed min_edge = 0.6 ,node_label_size = 1.2 ,node_size = 30 ,edge_size = 3 ,node_frame_col = "grey20" ,node_fill_col = "grey80" ,vertex.frame.width = 3
Visualizing a simplified representation of the circadian GCN
Annotate the network
Obtain a list of genes in each module
Show code
# module_genes <- purrr::map( # sample.names, # ~ mods[[.x]][["module_genes"]] # ) |> # setNames(sample.names) <- mods[["module_genes" ]]
Identify rhythmic modules
Show code
<- purrr:: map (~ load_rhy_genes (sample = .x|> setNames (sample.names)# l_module_genes <- purrr::map( # sample.names, # ~ module_genes[[.x]] |> # arrange(module_identity) |> # group_split(module_identity) |> # purrr::map( # ~ .x |> pull(gene_name) # ) |> # setNames(unique(module_genes[[.x]]$module_identity)) # ) |> # setNames(sample.names) <- module_genes |> arrange (module_identity) |> group_split (module_identity) |> :: map (~ .x |> pull (gene_name)|> setNames (unique (module_genes$ module_identity))###- ### - ### - ### - ### - ### - ### - ### - # Set your p-value of choice # col_pval = "BH.Q" = "default.pvalue" # col_pval = "raw.pvalue" ###- ### - ### - ### - ### - ### - ### - ### - <- purrr:: map (~ db_rhy[[.x]] |> :: map (~ .x |> filter (if_all (all_of (col_pval),~ .x < 0.05 |> filter (%in% unlist (l_module_genes)|> pull (1 ) |> unique ()|> :: compact ()|> setNames (sample.names)
Modules vs. Rhythmic genes
Show code
# LIST ONE - WGCNA modules <- l_module_genessapply (list1, length) |> print ()
C1 C2 C3 C4 C5 C6 C7 C8 C9 C10 C11 C12 C13 C14
154 163 3195 105 748 146 54 796 2462 440 194 61 50 325
Show code
<- purrr:: map (function (x) {cat ("Tissue:" , x, " \n " )## LIST TWO - rhythmic genes <- l_rhy_genes[[x]]sapply (list2, length) |> print ()## CHECK FOR OVERLAP # define size of genome = length (unique (c (unlist (list1), unlist (list2))))# make a GOM object .1 v2 <- GeneOverlap:: newGOM (genome.size = size:: drawHeatmap (.1 v2,adj.p = TRUE ,cutoff = 0.05 ,what = "odds.ratio" ,# what="Jaccard", log.scale = T,note.col = "black" ,grid.col = "Oranges" .1 v2
Tissue: mmus-brain_stem
ARS empJTK GeneCycle JTK meta2d RAIN
1010 361 368 198 783 687
Tissue: mmus-cerebellum
ARS empJTK GeneCycle JTK meta2d RAIN
2663 701 641 241 1952 593
Tissue: mmus-hypothalamus
ARS empJTK GeneCycle JTK meta2d RAIN
2042 559 481 223 1294 609
Tissue: mmus-heart
ARS empJTK GeneCycle JTK meta2d RAIN
2550 987 785 532 2251 1316
